Translating intracellular calcium signaling into models

The rich experimental data on intracellular calcium has put theoreticians in an ideal position to derive models of intracellular calcium signaling. Over the last 25 years, a large number of modeling frameworks have been suggested. Here, I will review some of the milestones of intracellular calcium modeling with a special emphasis on calcium-induced calcium release (CICR) through inositol-1,4,5-trisphosphate and ryanodine receptors. I will highlight key features of CICR and how they are represented in models as well as the challenges that theoreticians face when translating our current understanding of calcium signals into equations. The selected examples demonstrate that a successful model provides mechanistic insights into the molecular machinery of the Ca2+ signaling toolbox and determines the contribution of local Ca2+ release to global Ca2+ patterns, which at the moment cannot be resolved experimentally

Time to blip – stochastic simulation of single channel opening

The stochastic dynamics of the inositol-1,4,5-trisphosphate (IP3) receptor (IP3R) is key to understanding a wide range of observed calcium (Ca2+) signals (Falcke 2004). The stochastic nature results from the constant binding and unbinding of Ca2+ and IP3 to and from their respective binding sites and is especially important in the initiation of a Ca2+ puff, i.e. the release of Ca2+ through a cluster of IP3Rs. Once the first IP3R opens, the Ca2+ concentration rises significantly around the ion channel and hence increases the open probability for neighboring IP3Rs. In turn this may trigger the activation of further receptors giving rise to a Ca2+ puff (Thul et al. 2009; Thurley et al. 2012). In this protocol, we determine the time that it takes for a single IP3R to open from rest. We explicitly take into account the tetrameric structure of the IP3R and the fact that multiple subunits need to be active before the channel opens (Bezprozvanny et al. 1991; Watras et al. 1991). We develop code for a stochastic simulation of the IP3R and simulate it using the software package Matlab (Attaway 2011). This protocol demonstrates the basic form of a stochastic simulation algorithm and may serve as a starting point to investigate more complex gating dynamics

Oscillations in a point models of the intracellular Ca2+ concentration

Oscillations in the intracellular calcium (Ca2+) concentration form one of the main pathways by which cells translate external stimuli into physiological responses (Thul et al. 2008; Dupont et al. 2011; Parekh 2011). The mechanisms that underlie the generation of Ca2+ oscillations are still actively debated in the modeling community, but there is growing evidence that Ca2+ oscillations result from the spatio-temporal summation of subcellular Ca2+ release events (Thurley et al. 2012). Nevertheless, one prominent modeling approach to intracellular Ca2+ oscillations is the use of ordinary differential equations (ODEs), which treat the intracellular Ca2+ concentration as spatially homogenous. Although ODEs cannot account for the interaction of Ca2+ microdomains to form cell-wide Ca2+ patterns, modelers still choose ODEs since (a) the study of ODEs is computationally cheap, and a large body of techniques is available to investigate ODEs in great detail, or (b) there might not be sufficient experimental data to develop a spatially extended model. Irrespective of the reason, analyzing ODEs is a key instrument in the toolbox of modelers. In this protocol, we look at a wellknown model for Ca2+ oscillations (De Young and Keizer 1992; Li and Rinzel 1994). The main emphasis of this protocol is the use of the open source software package XPPaut to numerically study ODEs (Ermentrout 2002). The knowledge gained here can be directly transferred to other ODE systems and therefore may serve as a template for future studies. For a general background on analysing ODEs in the context of Mathematical Cell Physiology, I refer the reader to (Keener and Sneyd 2001; Fall et al. 2002; Britton 2002; Murray 2013)

Evolution of moments and correlations in non-renewal escape-time processes

The theoretical description of non-renewal stochastic systems is a challenge. Analytical results are often not available or can only be obtained under strong conditions, limiting their applicability. Also, numerical results have mostly been obtained by ad-hoc Monte--Carlo simulations, which are usually computationally expensive when a high degree of accuracy is needed. To gain quantitative insight into these systems under general conditions, we here introduce a numerical iterated first-passage time approach based on solving the time-dependent Fokker--Planck equation (FPE) to describe the statistics of non-renewal stochastic systems. We illustrate the approach using spike-triggered neuronal adaptation in the leaky and perfect integrate-and-fire model, respectively. The transition to stationarity of first-passage time moments and their sequential correlations occur on a non-trivial timescale that depends on all system parameters. Surprisingly this is so for both single exponential and scale-free power-law adaptation. The method works beyond the small noise and timescale separation approximations. It shows excellent agreement with direct Monte Carlo simulations, which allows for the computation of transient and stationary distributions. We compare different methods to compute the evolution of the moments and serial correlation coefficients (SCC), and discuss the challenge of reliably computing the SCC which we find to be very sensitive to numerical inaccuracies for both the leaky and perfect integrate-and-fire models. In conclusion, our methods provide a general picture of non-renewal dynamics in a wide range of stochastic systems exhibiting short and long-range correlations

First passage times in integrate-and-fire neurons with stochastic thresholds

We consider a leaky integrate--and--fire neuron with deterministic subthreshold dynamics and a firing threshold that evolves as an Ornstein-Uhlenbeck process. The formulation of this minimal model is motivated by the experimentally observed widespread variation of neural firing thresholds. We show numerically that the mean first passage time can depend non-monotonically on the noise amplitude. For sufficiently large values of the correlation time of the stochastic threshold the mean first passage time is maximal for non-vanishing noise. We provide an explanation for this effect by analytically transforming the original model into a first passage time problem for Brownian motion. This transformation also allows for a perturbative calculation of the first passage time histograms. In turn this provides quantitative insights into the mechanisms that lead to the non-monotonic behaviour of the mean first passage time. The perturbation expansion is in excellent agreement with direct numerical simulations. The approach developed here can be applied to any deterministic subthreshold dynamics and any Gauss-Markov processes for the firing threshold. This opens up the possibility to incorporate biophysically detailed components into the subthreshold dynamics, rendering our approach a powerful framework that sits between traditional integrate-and-fire models and complex mechanistic descriptions of neural dynamics

Networks of piecewise linear neural mass models

Neural mass models are ubiquitous in large scale brain modelling. At the node level they are written in terms of a set of ordinary differential equations with a nonlinearity that is typically a sigmoidal shape. Using structural data from brain atlases they may be connected into a network to investigate the emergence of functional dynamic states, such as synchrony. With the simple restriction of the classic sigmoidal nonlinearity to a piecewise linear caricature we show that the famous Wilson-Cowan neural mass model can be explicitly analysed at both the node and network level. The construction of periodic orbits at the node level is achieved by patching together matrix exponential solutions, and stability is determined using Floquet theory. For networks with interactions described by circulant matrices, we show that the stability of the synchronous state can be determined in terms of a low-dimensional Floquet problem parameterised by the eigenvalues of the interaction matrix. Moreover, this network Floquet problem is readily solved using linear algebra, to predict the onset of spatio-temporal network patterns arising from a synchronous instability. We further consider the case of a discontinuous choice for the node nonlinearity, namely the replacement of the sigmoid by a Heaviside nonlinearity. This gives rise to a continuous-time switching network. At the node level this allows for the existence of unstable sliding periodic orbits, which we explicitly construct. The stability of a periodic orbit is now treated with a modification of Floquet theory to treat the evolution of small perturbations through switching manifolds via the use of saltation matrices. At the network level the stability analysis of the synchronous state is considerably more challenging. Here we report on the use of ideas originally developed for the study of Glass networks to treat the stability of periodic network states in neural mass models with discontinuous interactions

Three-dimensional spatio-temporal modelling of store operated Ca2+ entry: insights into ER refilling and the spatial signature of Ca2+ signals

The spatial organisation of Orai channels and SERCA pumps within ER-PM junctions is important for enhancing the versatility and specificity of subcellular Ca2+ signals generated during store operated Ca2+ entry (SOCE). In this paper we present a novel three dimensional spatio-temporal model describing Ca2+ dynamics in the ER-PM junction and sub-PM ER during SOCE. We investigate the role of Orai channel and SERCA pump location to provide insights into how these components shape the Ca2+ signals generated and affect ER refilling. We find that the organisation of Orai channels within the ER-PM junction controls the amplitude and shape of the Ca2+ profile but does not enhance ER refilling. The model shows that ER refilling is only weakly affected by the location of SERCA2b pumps within the ER-PM junction and that the placement of SERCA2a pumps within the ER-PM junction has much greater control over ER refilling

A statistical view on calcium oscillations

Transient rises and falls of the intracellular calcium concentration have been observed in numerous cell types and under a plethora of conditions. There is now a growing body of evidence that these whole-cell calcium oscillations are stochastic, which poses a significant challenge for modelling. In this review, we take a closer look at recently developed statistical approaches to calcium oscillations. These models describe the timing of whole-cell calcium spikes, yet their parametrisations reflect subcellular processes. We show how non-stationary calcium spike sequences, which e.g. occur during slow depletion of intracellular calcium stores or in the presence of time-dependent stimulation, can be analysed with the help of so-called intensity functions. By utilising Bayesian concepts, we demonstrate how values of key parameters of the statistical model can be inferred from single cell calcium spike sequences and illustrate what information whole-cell statistical models can provide about the subcellular mechanistic processes that drive calcium oscillations. In particular, we find that the interspike interval distribution of HEK293 cells under constant stimulation is captured by a Gamma distribution

Release currents of IP₃ receptor channel clusters and concentration profiles

We simulate currents and concentration profiles generated by Ca2+ release from the endoplasmic reticulum (ER) to the cytosol through IP3 receptor channel clusters. Clusters are described as conducting pores in the lumenal membrane with a diameter from 6 nm to 36 nm. The endoplasmic reticulum is modeled as a disc with a radius of 1–12 mm and an inner height of 28 nm. We adapt the dependence of the currents on the trans Ca2+ concentration (intralumenal) measured in lipid bilayer experiments to the cellular geometry. Simulated currents are compared with signal mass measurements in Xenopus oocytes. We find that release currents depend linearly on the concentration of free Ca2+ in the lumen. The release current is approximately proportional to the square root of the number of open channels in a cluster. Cytosolic concentrations at the location of the cluster range from 25 μM to 170 μM. Concentration increase due to puffs in a distance of a few micrometers from the puff site is found to be in the nanomolar range. Release currents decay biexponentially with timescales of < 1 s and a few seconds. Concentration profiles decay with timescales of 0.125–0.250 s upon termination of release

Understanding cardiac alternans: a piecewise linear modeling framework

Cardiac alternans is a beat-to-beat alternation in action potential duration (APD) and intracellular calcium (Ca²⁺) cycling seen in cardiac myocytes under rapid pacing that is believed to be a precursor to fibrillation. The cellular mechanisms of these rhythms and the coupling between cellular Ca²⁺ and voltage dynamics have been extensively studied leading to the development of a class of physiologically detailed models. These have been shown numerically to reproduce many of the features of myocyte response to pacing, including alternans, and have been analyzed mathematically using various approximation techniques that allow for the formulation of a low dimensional map to describe the evolution of APDs. The seminal work by Shiferaw and Karma is of particular interest in this regard [Shiferaw, Y. and Karma, A., "Turing instability mediated by voltage and calcium diffusion in paced cardiac cells," Proc. Natl. Acad. Sci. U.S.A. 103, 5670-5675 (2006)]. Here, we establish that the key dynamical behaviors of the Shiferaw-Karma model are arranged around a set of switches. These are shown to be the main elements for organizing the nonlinear behavior of the model. Exploiting this observation, we show that a piecewise linear caricature of the Shiferaw-Karma model, with a set of appropriate switching manifolds, can be constructed that preserves the physiological interpretation of the original model while being amenable to a systematic mathematical analysis. In illustration of this point, we formulate the dynamics of Ca²⁺ cycling (in response to pacing) and compute the properties of periodic orbits in terms of a stroboscopic map that can be constructed without approximation. Using this, we show that alternans emerge via a period-doubling instability and track this bifurcation in terms of physiologically important parameters. We also show that when coupled to a spatially extended model for Ca²⁺ transport, the model supports spatially varying patterns of alternans. We analyze the onset of this instability with a generalization of the master stability approach to accommodate the nonsmooth nature of our system